Setup
knitr::opts_chunk$set(cache = TRUE)
library(tidyverse)
library(leaflet)
library(sf)
library(dbscan)
library(spatstat)
library(splancs)
library(htmlwidgets)
library(webshot)
data_dir <- "../clean_data/"
maps_dir <- "../maps/"
EDA & Data Wrangling
starbucks <- read_csv(paste0(data_dir, "starbucks_shanghai.csv")) |>
st_as_sf(coords = c("longitude", "latitude"), crs = 4326) |>
select(-city)
## Rows: 754 Columns: 3
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): city
## dbl (2): latitude, longitude
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
starbucks
## Simple feature collection with 754 features and 0 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 121.0167 ymin: 30.73263 xmax: 121.9256 ymax: 31.62088
## Geodetic CRS: WGS 84
## # A tibble: 754 × 1
## geometry
## <POINT [°]>
## 1 (121.3965 31.62088)
## 2 (121.2479 31.39764)
## 3 (121.2374 31.39053)
## 4 (121.3506 31.40225)
## 5 (121.2456 31.38443)
## 6 (121.2562 31.38369)
## 7 (121.2788 31.38513)
## 8 (121.3364 31.39046)
## 9 (121.3583 31.38815)
## 10 (121.2587 31.35715)
## # ℹ 744 more rows
shanghai <- st_read(paste0(maps_dir, "2019年县级.shp")) |>
filter(NAME_1 == "Shanghai") |>
select(ENG_NAME) |>
rename(District = ENG_NAME) |>
mutate(District = case_when(
District == "Pudongxin" ~ "Pudong New Area",
District == "Jingan" ~ "Jing'an",
TRUE ~ District
))
## Reading layer `2019年县级' from data source
## `/Users/bulldf/University of Toronto/2024 Fall/STA465/project/maps/2019年县级.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 2856 features and 29 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 73.50114 ymin: 3.837902 xmax: 135.0885 ymax: 53.5609
## Geodetic CRS: WGS 84
shanghai
## Simple feature collection with 16 features and 1 field
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 120.8523 ymin: 30.6917 xmax: 122.1129 ymax: 31.87255
## Geodetic CRS: WGS 84
## First 10 features:
## District geometry
## 1 Chongming MULTIPOLYGON (((121.4272 31...
## 2 Fengxian MULTIPOLYGON (((121.5661 31...
## 3 Hongkou MULTIPOLYGON (((121.4808 31...
## 4 Huangpu MULTIPOLYGON (((121.4904 31...
## 5 Jiading MULTIPOLYGON (((121.3318 31...
## 6 Jinshan MULTIPOLYGON (((121.0238 30...
## 7 Jing'an MULTIPOLYGON (((121.4638 31...
## 8 Minhang MULTIPOLYGON (((121.3336 31...
## 9 Pudong New Area MULTIPOLYGON (((121.8514 31...
## 10 Qingpu MULTIPOLYGON (((121.1492 31...
demographics <- read_csv(paste0(data_dir, "shanghai_2019_en.csv"))
## Rows: 16 Columns: 4
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): District
## dbl (3): Area, Population, Density
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
demographics
## # A tibble: 16 × 4
## District Area Population Density
## <chr> <dbl> <dbl> <dbl>
## 1 Pudong New Area 1210. 557. 4599
## 2 Huangpu 20.5 65.1 31808
## 3 Xuhui 54.8 109. 19989
## 4 Changning 38.3 69.4 18110
## 5 Jing'an 36.9 106. 28680
## 6 Putuo 54.8 128. 23268
## 7 Hongkou 23.5 79.4 33816
## 8 Yangpu 60.7 130. 21487
## 9 Minhang 371. 255. 6876
## 10 Baoshan 271. 204. 7544
## 11 Jiading 464. 160. 3438
## 12 Jinshan 586. 80.7 1377
## 13 Songjiang 606. 177. 2926
## 14 Qingpu 670. 123. 1840
## 15 Fengxian 687. 116. 1684
## 16 Chongming 1185. 68.4 577
leaflet() |>
addProviderTiles("CartoDB.Positron") |>
addCircles(
data = starbucks,
color = "darkgreen"
) |>
addPolygons(
data = shanghai,
color = "#4B4B4B",
weight = 1.5,
opacity = 1,
label = ~District
)
wgs_coords <- st_coordinates(starbucks)
utm_coords <- starbucks |> st_transform(crs = 32651) |> st_coordinates()
df <- starbucks |>
mutate(
lon = wgs_coords[, 1],
lat = wgs_coords[, 2],
x = utm_coords[, 1],
y = utm_coords[, 2]
) |>
st_drop_geometry()
df
## # A tibble: 754 × 4
## lon lat x y
## * <dbl> <dbl> <dbl> <dbl>
## 1 121. 31.6 347914. 3499530.
## 2 121. 31.4 333421. 3474999.
## 3 121. 31.4 332409. 3474227.
## 4 121. 31.4 343190. 3475359.
## 5 121. 31.4 333176. 3473538.
## 6 121. 31.4 334182. 3473440.
## 7 121. 31.4 336336. 3473566.
## 8 121. 31.4 341819. 3474072.
## 9 121. 31.4 343900. 3473785.
## 10 121. 31.4 334376. 3470495.
## # ℹ 744 more rows
window <- convexhull.xy(df$x, df$y)
window
## window: polygonal boundary
## enclosing rectangle: [310427.1, 397312.9] x [3401137, 3499530] units
starbucks_ppp <- ppp(df$x, df$y, window = window)
starbucks_ppp
## Planar point pattern: 754 points
## window: polygonal boundary
## enclosing rectangle: [310427.1, 397312.9] x [3401137, 3499530] units
starbucks_pts <- as.points(df$x, df$y)
plot(starbucks_ppp, main = "Point Pattern of Starbucks Stores in Shanghai")
Homogeneous PPP
hppp <- rpoispp(lambda = nrow(df) / area(window), win = window)
hppp
## Planar point pattern: 736 points
## window: polygonal boundary
## enclosing rectangle: [310427.1, 397312.9] x [3401137, 3499530] units
plot(hppp)
Testing for CSR
Ripley’s K
min_x <- min(df$x)
max_x <- max(df$x)
min_y <- min(df$y)
max_y <- max(df$y)
poly <- as.points(c(min_x, max_x, max_x, min_x), c(min_y, min_y, max_y, max_y))
starbucks_seq <- seq(0, 80000, 2000)
khat <- khat(starbucks_pts, poly, starbucks_seq)
khat
## [1] 0 242653037 790222405 1510179002 2284491196 3072215026
## [7] 3797145392 4442974302 5023099746 5524568337 5952076682 6312154056
## [13] 6623368200 6893070622 7142798921 7378406327 7608843623 7840342986
## [19] 8048979070 8217006776 8363620828 8489521559 8598801337 8688037900
## [25] 8780845469 8875149534 8967234186 9052615422 9174949555 9309131136
## [31] 9409884731 9489957803 9557019064 9609670451 9651253736 9678783948
## [37] 9694863508 9709981396 9718607706 9723013929 9729432091
ul_khat <- Kenv.csr(length(df$x), poly, nsim = 99, starbucks_seq)
## Doing simulation 1
## Doing simulation 2
## Doing simulation 3
## Doing simulation 4
## Doing simulation 5
## Doing simulation 6
## Doing simulation 7
## Doing simulation 8
## Doing simulation 9
## Doing simulation 10
## Doing simulation 11
## Doing simulation 12
## Doing simulation 13
## Doing simulation 14
## Doing simulation 15
## Doing simulation 16
## Doing simulation 17
## Doing simulation 18
## Doing simulation 19
## Doing simulation 20
## Doing simulation 21
## Doing simulation 22
## Doing simulation 23
## Doing simulation 24
## Doing simulation 25
## Doing simulation 26
## Doing simulation 27
## Doing simulation 28
## Doing simulation 29
## Doing simulation 30
## Doing simulation 31
## Doing simulation 32
## Doing simulation 33
## Doing simulation 34
## Doing simulation 35
## Doing simulation 36
## Doing simulation 37
## Doing simulation 38
## Doing simulation 39
## Doing simulation 40
## Doing simulation 41
## Doing simulation 42
## Doing simulation 43
## Doing simulation 44
## Doing simulation 45
## Doing simulation 46
## Doing simulation 47
## Doing simulation 48
## Doing simulation 49
## Doing simulation 50
## Doing simulation 51
## Doing simulation 52
## Doing simulation 53
## Doing simulation 54
## Doing simulation 55
## Doing simulation 56
## Doing simulation 57
## Doing simulation 58
## Doing simulation 59
## Doing simulation 60
## Doing simulation 61
## Doing simulation 62
## Doing simulation 63
## Doing simulation 64
## Doing simulation 65
## Doing simulation 66
## Doing simulation 67
## Doing simulation 68
## Doing simulation 69
## Doing simulation 70
## Doing simulation 71
## Doing simulation 72
## Doing simulation 73
## Doing simulation 74
## Doing simulation 75
## Doing simulation 76
## Doing simulation 77
## Doing simulation 78
## Doing simulation 79
## Doing simulation 80
## Doing simulation 81
## Doing simulation 82
## Doing simulation 83
## Doing simulation 84
## Doing simulation 85
## Doing simulation 86
## Doing simulation 87
## Doing simulation 88
## Doing simulation 89
## Doing simulation 90
## Doing simulation 91
## Doing simulation 92
## Doing simulation 93
## Doing simulation 94
## Doing simulation 95
## Doing simulation 96
## Doing simulation 97
## Doing simulation 98
## Doing simulation 99
plot(
starbucks_seq,
khat - pi * starbucks_seq^2,
type = "l",
xlab = "Distance",
ylab = "Estimated K - pi * h^2",
main = "Ripley's K"
)
# plot upper bound
lines(starbucks_seq, ul_khat$upper - pi * starbucks_seq^2, lty = 2)
# plot lower bound
lines(starbucks_seq, ul_khat$lower - pi * starbucks_seq^2, lty = 2)
l <- function(k, h) {
sqrt(k / pi) - h
}
plot(
starbucks_seq,
l(khat, starbucks_seq),
type = "l",
xlab = "Distance",
ylab = "Estimated L"
)
# plot upper bound of Lhat
lines(starbucks_seq, l(ul_khat$upper, starbucks_seq), lty = 2)
# plot lower bound of Lhat
lines(starbucks_seq, l(ul_khat$lower, starbucks_seq), lty = 2)
K-S Test
ks_x <- cdf.test(starbucks_ppp, test = "ks", "x")
plot(ks_x)
ks_y <- cdf.test(starbucks_ppp, test = "ks", "y")
## Warning in ks.test.default(U, "punif", ...): ties should not be present for the
## Kolmogorov-Smirnov test
plot(ks_y)
G-function
plot(envelope(starbucks_ppp, Gest), main = "G-Function")
## Generating 99 simulations of CSR ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
## 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
## 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
## 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
## 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98,
## 99.
##
## Done.
Quadrat Counting
quadrat <- quadratcount(starbucks_ppp, nx = 10, ny = 10)
plot(
starbucks_pts,
cex = 0.5,
main = "Starbucks Quadrat Counts",
xlab = "x",
ylab = "y"
)
plot(quadrat, add = TRUE)
mean(quadrat)
## [1] 10.05333
quadrat.test(quadrat, alternative = "two.sided", method = "MonteCarlo")
##
## Conditional Monte Carlo test of CSR using quadrat counts
## Test statistic: Pearson X2 statistic
##
## data:
## X2 = 4821.2, p-value = 0.001
## alternative hypothesis: two.sided
##
## Quadrats: 75 tiles (irregular windows)
quadrat.test(quadrat, alternative = "regular", method = "MonteCarlo")
##
## Conditional Monte Carlo test of CSR using quadrat counts
## Test statistic: Pearson X2 statistic
##
## data:
## X2 = 4821.2, p-value = 1
## alternative hypothesis: regular
##
## Quadrats: 75 tiles (irregular windows)
quadrat.test(quadrat, alternative = "clustered", method = "MonteCarlo")
##
## Conditional Monte Carlo test of CSR using quadrat counts
## Test statistic: Pearson X2 statistic
##
## data:
## X2 = 4821.2, p-value = 5e-04
## alternative hypothesis: clustered
##
## Quadrats: 75 tiles (irregular windows)
Kernel Density
optim_bw <- bw.diggle(starbucks_ppp, edge = "border")
optim_bw
## sigma
## 340.0615
ds1 <- density.ppp(starbucks_ppp, optim_bw)
plot(ds1, main = "Starbucks Density for Optimal Bandwidth")
plot(starbucks_ppp, add = TRUE, cols = "white", pch = 16, size = 0.5)
ds2 <- density.ppp(starbucks_ppp, sigma = 3000)
plot(ds2, main = "Starbucks Density for Bandwidth = 3000")
plot(starbucks_ppp, add = TRUE, cols = "white", pch = 16, size = 0.5)
plot(density(starbucks_ppp), main = "Starbucks Density")
plot(starbucks_ppp, add = TRUE, cols = "white", pch = 16, size = 0.5)
Dirichlet Tesselation
plot(dirichlet(starbucks_ppp), main = "Starbucks Dirichlet Tesselation")
DBSCAN
dbscan1 <- dbscan(starbucks_pts, eps = 500, minPts = 5)
dbscan1
## DBSCAN clustering for 754 objects.
## Parameters: eps = 500, minPts = 5
## Using euclidean distances and borderpoints = TRUE
## The clustering contains 16 cluster(s) and 539 noise points.
##
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
## 539 9 8 5 45 8 16 44 15 6 6 2 5 26 4 5 11
##
## Available fields: cluster, eps, minPts, metric, borderPoints
hullplot(starbucks_pts, dbscan1$cluster, main = "DBSCAN Clusters")
## Warning in hullplot(starbucks_pts, dbscan1$cluster, main = "DBSCAN Clusters"):
## Not enough colors. Some colors will be reused.
dbscan2 <- dbscan(starbucks_pts, eps = 2000, minPts = 5)
dbscan2
## DBSCAN clustering for 754 objects.
## Parameters: eps = 2000, minPts = 5
## Using euclidean distances and borderpoints = TRUE
## The clustering contains 9 cluster(s) and 134 noise points.
##
## 0 1 2 3 4 5 6 7 8 9
## 134 5 536 20 6 10 15 8 5 15
##
## Available fields: cluster, eps, minPts, metric, borderPoints
hullplot(starbucks_pts, dbscan2$cluster, main = "DBSCAN Clusters")
## Warning in hullplot(starbucks_pts, dbscan2$cluster, main = "DBSCAN Clusters"):
## Not enough colors. Some colors will be reused.
dbscan3 <- dbscan(starbucks_pts, eps = 5000, minPts = 5)
dbscan3
## DBSCAN clustering for 754 objects.
## Parameters: eps = 5000, minPts = 5
## Using euclidean distances and borderpoints = TRUE
## The clustering contains 7 cluster(s) and 23 noise points.
##
## 0 1 2 3 4 5 6 7
## 23 674 6 10 15 5 6 15
##
## Available fields: cluster, eps, minPts, metric, borderPoints
hullplot(starbucks_pts, dbscan3$cluster, main = "DBSCAN Clusters")
dbscan4 <- hdbscan(starbucks_pts, minPts = 5)
dbscan4
## HDBSCAN clustering for 754 objects.
## Parameters: minPts = 5
## The clustering contains 43 cluster(s) and 320 noise points.
##
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
## 320 5 6 6 10 15 14 6 6 10 5 26 8 12 15 10 7 10 6 18
## 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39
## 7 10 6 6 5 7 19 8 32 10 10 8 5 7 10 20 15 5 10 5
## 40 41 42 43
## 13 7 6 8
##
## Available fields: cluster, minPts, coredist, cluster_scores,
## membership_prob, outlier_scores, hc
hullplot(starbucks_pts, dbscan4$cluster, main = "HDBSCAN Clusters")
## Warning in hullplot(starbucks_pts, dbscan4$cluster, main = "HDBSCAN Clusters"):
## Not enough colors. Some colors will be reused.
plot(dbscan4, show_flat = TRUE)
Inhomogeneous Poisson Process
model1 <- ppm(starbucks_ppp ~ 1, Poisson())
summary(model1)
## Point process model
## Fitted to data: starbucks_ppp
## Fitting method: maximum likelihood
## Model was fitted analytically
## Call:
## ppm.formula(Q = starbucks_ppp ~ 1, interaction = Poisson())
## Edge correction: "border"
## [border correction distance r = 0 ]
## --------------------------------------------------------------------------------
## Quadrature scheme (Berman-Turner) = data + dummy + weights
##
## Data pattern:
## Planar point pattern: 754 points
## Average intensity 1.54e-07 points per square unit
## Window: polygonal boundary
## single connected closed polygon with 10 vertices
## enclosing rectangle: [310427.1, 397312.9] x [3401137, 3499530] units
## (86890 x 98390 units)
## Window area = 4894750000 square units
## Fraction of frame area: 0.573
##
## Dummy quadrature points:
## 64 x 64 grid of dummy points, plus 4 corner points
## dummy spacing: 1357.589 x 1537.401 units
##
## Original dummy parameters: =
## Planar point pattern: 2410 points
## Average intensity 4.92e-07 points per square unit
## Window: polygonal boundary
## single connected closed polygon with 10 vertices
## enclosing rectangle: [310427.1, 397312.9] x [3401137, 3499530] units
## (86890 x 98390 units)
## Window area = 4894750000 square units
## Fraction of frame area: 0.573
## Quadrature weights:
## (counting weights based on 64 x 64 array of rectangular tiles)
## All weights:
## range: [110000, 2090000] total: 4.89e+09
## Weights on data points:
## range: [110000, 1040000] total: 3.95e+08
## Weights on dummy points:
## range: [110000, 2090000] total: 4.49e+09
## --------------------------------------------------------------------------------
## FITTED :
##
## Stationary Poisson process
##
## ---- Intensity: ----
##
##
## Uniform intensity:
## [1] 1.540426e-07
##
## Estimate S.E. CI95.lo CI95.hi Ztest Zval
## log(lambda) -15.68604 0.03641785 -15.75741 -15.61466 *** -430.7238
##
## ----------- gory details -----
##
## Fitted regular parameters (theta):
## log(lambda)
## -15.68604
##
## Fitted exp(theta):
## log(lambda)
## 1.540426e-07
AIC(model1)
## [1] 25164.54
model2 <- ppm(starbucks_ppp, ~ log(x) + log(y), Poisson())
summary(model2)
## Point process model
## Fitted to data: starbucks_ppp
## Fitting method: maximum likelihood (Berman-Turner approximation)
## Model was fitted using glm()
## Algorithm converged
## Call:
## ppm.ppp(Q = starbucks_ppp, trend = ~log(x) + log(y), interaction = Poisson())
## Edge correction: "border"
## [border correction distance r = 0 ]
## --------------------------------------------------------------------------------
## Quadrature scheme (Berman-Turner) = data + dummy + weights
##
## Data pattern:
## Planar point pattern: 754 points
## Average intensity 1.54e-07 points per square unit
## Window: polygonal boundary
## single connected closed polygon with 10 vertices
## enclosing rectangle: [310427.1, 397312.9] x [3401137, 3499530] units
## (86890 x 98390 units)
## Window area = 4894750000 square units
## Fraction of frame area: 0.573
##
## Dummy quadrature points:
## 64 x 64 grid of dummy points, plus 4 corner points
## dummy spacing: 1357.589 x 1537.401 units
##
## Original dummy parameters: =
## Planar point pattern: 2410 points
## Average intensity 4.92e-07 points per square unit
## Window: polygonal boundary
## single connected closed polygon with 10 vertices
## enclosing rectangle: [310427.1, 397312.9] x [3401137, 3499530] units
## (86890 x 98390 units)
## Window area = 4894750000 square units
## Fraction of frame area: 0.573
## Quadrature weights:
## (counting weights based on 64 x 64 array of rectangular tiles)
## All weights:
## range: [110000, 2090000] total: 4.89e+09
## Weights on data points:
## range: [110000, 1040000] total: 3.95e+08
## Weights on dummy points:
## range: [110000, 2090000] total: 4.49e+09
## --------------------------------------------------------------------------------
## FITTED :
##
## Nonstationary Poisson process
##
## ---- Intensity: ----
##
## Log intensity: ~log(x) + log(y)
##
## Fitted trend coefficients:
## (Intercept) log(x) log(y)
## -1400.234361 3.253521 89.217040
##
## Estimate S.E. CI95.lo CI95.hi Ztest Zval
## (Intercept) -1400.234361 87.7026941 -1572.128483 -1228.340239 *** -15.965694
## log(x) 3.253521 0.7566339 1.770546 4.736497 *** 4.299995
## log(y) 89.217040 5.6957551 78.053565 100.380515 *** 15.663777
##
## ----------- gory details -----
##
## Fitted regular parameters (theta):
## (Intercept) log(x) log(y)
## -1400.234361 3.253521 89.217040
##
## Fitted exp(theta):
## (Intercept) log(x) log(y)
## 0.000000e+00 2.588132e+01 5.577866e+38
AIC(model2)
## [1] 24909.49
plot(model2)
simulation1 <- simulate(model1, nsim = 3)
## Generating 3 simulated patterns ...1, 2,
## 3.
simulation1[[4]] <- starbucks_ppp
plot(simulation1, main = "")
simulation2 <- simulate(model2, nsim = 3)
## Generating 3 simulated patterns ...1, 2,
## 3.
simulation2[[4]] <- starbucks_ppp
plot(simulation2, main = "")
Cluster Process
model3 <- kppm(starbucks_ppp, clusters = "Thomas")
summary(model3)
## Stationary cluster point process model
## Fitted to point pattern dataset 'starbucks_ppp'
## Fitted by minimum contrast
## Summary statistic: K-function
## Minimum contrast fit (object of class "minconfit")
## Model: Thomas process
## Fitted by matching theoretical K function to starbucks_ppp
##
## Internal parameters fitted by minimum contrast ($par):
## kappa sigma2
## 4.527876e-10 2.078070e+07
##
## Fitted cluster parameters:
## kappa scale
## 4.527876e-10 4.558585e+03
## Mean cluster size: 340.6065 points
##
## Converged successfully after 127 function evaluations
##
## Starting values of parameters:
## kappa sigma2
## 1.540426e-07 1.708917e+06
## Domain of integration: [ 0 , 21720 ]
## Exponents: p= 2, q= 0.25
##
## ----------- TREND -----
## Point process model
## Fitted to data: X
## Fitting method: maximum likelihood (Berman-Turner approximation)
## Model was fitted using glm()
## Algorithm converged
## Call:
## ppm.ppp(Q = X, trend = trend, rename.intercept = FALSE, covariates = covariates,
## covfunargs = covfunargs, use.gam = use.gam, forcefit = TRUE,
## improve.type = ppm.improve.type, improve.args = ppm.improve.args,
## nd = nd, eps = eps)
## Edge correction: "border"
## [border correction distance r = 0 ]
## --------------------------------------------------------------------------------
## Quadrature scheme (Berman-Turner) = data + dummy + weights
##
## Data pattern:
## Planar point pattern: 754 points
## Average intensity 1.54e-07 points per square unit
## Window: polygonal boundary
## single connected closed polygon with 10 vertices
## enclosing rectangle: [310427.1, 397312.9] x [3401137, 3499530] units
## (86890 x 98390 units)
## Window area = 4894750000 square units
## Fraction of frame area: 0.573
##
## Dummy quadrature points:
## 64 x 64 grid of dummy points, plus 4 corner points
## dummy spacing: 1357.589 x 1537.401 units
##
## Original dummy parameters: =
## Planar point pattern: 2410 points
## Average intensity 4.92e-07 points per square unit
## Window: polygonal boundary
## single connected closed polygon with 10 vertices
## enclosing rectangle: [310427.1, 397312.9] x [3401137, 3499530] units
## (86890 x 98390 units)
## Window area = 4894750000 square units
## Fraction of frame area: 0.573
## Quadrature weights:
## (counting weights based on 64 x 64 array of rectangular tiles)
## All weights:
## range: [110000, 2090000] total: 4.89e+09
## Weights on data points:
## range: [110000, 1040000] total: 3.95e+08
## Weights on dummy points:
## range: [110000, 2090000] total: 4.49e+09
## --------------------------------------------------------------------------------
## FITTED :
##
## Stationary Poisson process
##
## ---- Intensity: ----
##
##
## Uniform intensity:
## [1] 1.542224e-07
##
## Estimate S.E. CI95.lo CI95.hi Ztest Zval
## (Intercept) -15.68487 0.03641785 -15.75625 -15.61349 *** -430.6918
##
## ----------- gory details -----
##
## Fitted regular parameters (theta):
## (Intercept)
## -15.68487
##
## Fitted exp(theta):
## (Intercept)
## 1.542224e-07
##
## ----------- CLUSTER -----------
## Model: Thomas process
##
## Fitted cluster parameters:
## kappa scale
## 4.527876e-10 4.558585e+03
## Mean cluster size: 340.6065 points
##
## Final standard error and CI
## (allowing for correlation of cluster process):
## Estimate S.E. CI95.lo CI95.hi Ztest Zval
## (Intercept) -15.68487 0.6235393 -16.90698 -14.46276 *** -25.15458
##
## ----------- cluster strength indices ----------
## Sibling probability 0.8942624
## Count overdispersion index (on original window): 298.773
## Cluster strength: 8.457373
##
## Spatial persistence index (over window): 0
##
## Bound on distance from Poisson process (over window): 1
## = min (1, 1509.76, 5076500, 4819384, 2195.309)
##
## Bound on distance from MIXED Poisson process (over window): 1
##
## Intensity of parents of nonempty clusters: 4.527876e-10
## Mean number of offspring in a nonempty cluster: 340.6065
## Intensity of parents of clusters of more than one offspring point: 4.527876e-10
## Ratio of parents to parents-plus-offspring: 0.002927345 (where 1 = Poisson
## process)
## Probability that a typical point belongs to a nontrivial cluster: 1
model4 <- kppm(starbucks_ppp, clusters = "LGCP")
summary(model4)
## Stationary Cox point process model
## Fitted to point pattern dataset 'starbucks_ppp'
## Fitted by minimum contrast
## Summary statistic: K-function
## Minimum contrast fit (object of class "minconfit")
## Model: Log-Gaussian Cox process
## Covariance model: exponential
## Fitted by matching theoretical K function to starbucks_ppp
##
## Internal parameters fitted by minimum contrast ($par):
## sigma2 alpha
## 2.806734 11590.392857
##
## Fitted covariance parameters:
## var scale
## 2.806734 11590.392857
## Fitted mean of log of random intensity: -17.08824
##
## Converged successfully after 59 function evaluations
##
## Starting values of parameters:
## sigma2 alpha
## 1.000 1307.255
## Domain of integration: [ 0 , 21720 ]
## Exponents: p= 2, q= 0.25
##
## ----------- TREND -----
## Point process model
## Fitted to data: X
## Fitting method: maximum likelihood (Berman-Turner approximation)
## Model was fitted using glm()
## Algorithm converged
## Call:
## ppm.ppp(Q = X, trend = trend, rename.intercept = FALSE, covariates = covariates,
## covfunargs = covfunargs, use.gam = use.gam, forcefit = TRUE,
## improve.type = ppm.improve.type, improve.args = ppm.improve.args,
## nd = nd, eps = eps)
## Edge correction: "border"
## [border correction distance r = 0 ]
## --------------------------------------------------------------------------------
## Quadrature scheme (Berman-Turner) = data + dummy + weights
##
## Data pattern:
## Planar point pattern: 754 points
## Average intensity 1.54e-07 points per square unit
## Window: polygonal boundary
## single connected closed polygon with 10 vertices
## enclosing rectangle: [310427.1, 397312.9] x [3401137, 3499530] units
## (86890 x 98390 units)
## Window area = 4894750000 square units
## Fraction of frame area: 0.573
##
## Dummy quadrature points:
## 64 x 64 grid of dummy points, plus 4 corner points
## dummy spacing: 1357.589 x 1537.401 units
##
## Original dummy parameters: =
## Planar point pattern: 2410 points
## Average intensity 4.92e-07 points per square unit
## Window: polygonal boundary
## single connected closed polygon with 10 vertices
## enclosing rectangle: [310427.1, 397312.9] x [3401137, 3499530] units
## (86890 x 98390 units)
## Window area = 4894750000 square units
## Fraction of frame area: 0.573
## Quadrature weights:
## (counting weights based on 64 x 64 array of rectangular tiles)
## All weights:
## range: [110000, 2090000] total: 4.89e+09
## Weights on data points:
## range: [110000, 1040000] total: 3.95e+08
## Weights on dummy points:
## range: [110000, 2090000] total: 4.49e+09
## --------------------------------------------------------------------------------
## FITTED :
##
## Stationary Poisson process
##
## ---- Intensity: ----
##
##
## Uniform intensity:
## [1] 1.542224e-07
##
## Estimate S.E. CI95.lo CI95.hi Ztest Zval
## (Intercept) -15.68487 0.03641785 -15.75625 -15.61349 *** -430.6918
##
## ----------- gory details -----
##
## Fitted regular parameters (theta):
## (Intercept)
## -15.68487
##
## Fitted exp(theta):
## (Intercept)
## 1.542224e-07
##
## ----------- COX -----------
## Model: log-Gaussian Cox process
##
## Covariance model: exponential
## Fitted covariance parameters:
## var scale
## 2.806734 11590.392857
## Fitted mean of log of random intensity: -17.08824
##
## Final standard error and CI
## (allowing for correlation of Cox process):
## Estimate S.E. CI95.lo CI95.hi Ztest Zval
## (Intercept) -15.68487 0.7361998 -17.1278 -14.24195 *** -21.30518
simulation3 <- simulate(model3, nsim = 3)
simulation3[[4]] <- starbucks_ppp
plot(simulation3, main = "")
simulation4 <- simulate(model4, nsim = 3)
simulation4[[4]] <- starbucks_ppp
plot(simulation4, main = "")
